{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "hw10_anomaly_detection.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Iallen520/lhy_DL_Hw/blob/master/hw10_anomaly_detection.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YiVfKn-6tXz8",
        "colab_type": "text"
      },
      "source": [
        "# **Homework 10 - Anomaly Detection**\n",
        "\n",
        "若有任何問題，歡迎來信至助教信箱 ntu-ml-2020spring-ta@googlegroups.com"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6qLcLuhETSDH",
        "colab_type": "text"
      },
      "source": [
        "# Load and preprocess data"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cBo2oxu_WmZY",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!gdown --id '1_zT3JOpvXFGr7mkxs3XJDeGxTn_8pItq' --output train.npy \n",
        "!gdown --id '11Y_6JDjlhIY-M5-jW1rLRshDMqeKi9Kr' --output test.npy \n",
        "\n",
        "import numpy as np\n",
        "\n",
        "train = np.load('train.npy', allow_pickle=True)\n",
        "test = np.load('test.npy', allow_pickle=True)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LrpPP6x3XHo_",
        "colab_type": "text"
      },
      "source": [
        "# Task"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JbTcv510uhi6",
        "colab_type": "text"
      },
      "source": [
        "這份作業要執行的task是semi-supervised anomaly detection，也就是說trainingset是乾淨的，testing的時候才會混進outlier data(anomaly)。\n",
        "我們以某個簡單的image dataset（image加上他們的label（分類））作為示範，training data為原先trainingset中的某幾類，而testing data則是原先testingset的所有data，要偵測的anomaly為training data中未出現的類別。label的部分，1為outlier data而0為inlier data(相對於 outlier)。正確率以AUC計算。\n",
        "方法則列舉3種： KNN, PCA, Autoencoder"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QIHQqB7IXOXD",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "task = 'pca'"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q4GwvK4DXuW4",
        "colab_type": "text"
      },
      "source": [
        "# KNN"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qAIl7IETpbP2",
        "colab_type": "text"
      },
      "source": [
        "K-Nearest-Neighbor(KNN): 假設training data的label種類不多（e.g. < 20），然而因其為未知，可以猜測其為n，亦即假設training data有n群。先用K-means計算training data中的n個centroid，再用這n個centroid對training data分群。應該可以觀察到，inlier data與所分到群的centroid的距離應較outlier的此距離來得小。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "B7AB9wnaX1To",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from sklearn.cluster import MiniBatchKMeans\n",
        "from sklearn.metrics import f1_score, pairwise_distances, roc_auc_score\n",
        "from scipy.cluster.vq import vq, kmeans\n",
        "\n",
        "\n",
        "if task == 'knn':\n",
        "    x = train.reshape(len(train), -1)\n",
        "    y = test.reshape(len(test), -1)\n",
        "    scores = list()\n",
        "    for n in range(1, 10):\n",
        "      kmeans_x = MiniBatchKMeans(n_clusters=n, batch_size=100).fit(x)\n",
        "      y_cluster = kmeans_x.predict(y)\n",
        "      y_dist = np.sum(np.square(kmeans_x.cluster_centers_[y_cluster] - y), axis=1)\n",
        "\n",
        "      y_pred = y_dist\n",
        "      # score = f1_score(y_label, y_pred, average='micro')\n",
        "      # score = roc_auc_score(y_label, y_pred, average='micro')\n",
        "      # scores.append(score)\n",
        "    # print(np.max(scores), np.argmax(scores))\n",
        "    # print(scores)\n",
        "    # print('auc score: {}'.format(np.max(scores)))\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QLScetmgf4XV",
        "colab_type": "text"
      },
      "source": [
        "# PCA"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7ZKv6vNYMiWA",
        "colab_type": "text"
      },
      "source": [
        "PCA:首先計算training data的principle component，將testing data投影在這些component上，再將這些投影重建回原先space的向量。對重建的圖片和原圖計算MSE，inlier data的數值應該較outlier的數值為小。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EzgZGG7Pf6Qm",
        "colab_type": "code",
        "outputId": "5e38f33f-bca1-4e9a-8476-bdda926737a5",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "source": [
        "from sklearn.decomposition import PCA\n",
        "\n",
        "if task == 'pca':\n",
        "\n",
        "    x = train.reshape(len(train), -1)\n",
        "    y = test.reshape(len(test), -1)\n",
        "    pca = PCA(n_components=2).fit(x)\n",
        "\n",
        "    y_projected = pca.transform(y)\n",
        "    y_reconstructed = pca.inverse_transform(y_projected)  \n",
        "    dist = np.sqrt(np.sum(np.square(y_reconstructed - y).reshape(len(y), -1), axis=1))\n",
        "    \n",
        "    y_pred = dist\n",
        "    # score = roc_auc_score(y_label, y_pred, average='micro')\n",
        "    # score = f1_score(y_label, y_pred, average='micro')\n",
        "    # print('auc score: {}'.format(score))"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "auc score: 0.568625\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zR9zC0_Df-CR",
        "colab_type": "text"
      },
      "source": [
        "# Autoencoder"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1EbfwRREhA7c",
        "colab_type": "text"
      },
      "source": [
        "# Models & loss"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xUUIwFm02qEL",
        "colab_type": "text"
      },
      "source": [
        "課程影片參見：https://www.youtube.com/watch?v=6W8FqUGYyDo&list=PLJV_el3uVTsOK_ZK5L0Iv_EQoL1JefRL4&index=8"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lESPM7C9svYa",
        "colab_type": "text"
      },
      "source": [
        "fcn_autoencoder and vae are from https://github.com/L1aoXingyu/pytorch-beginner"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GwG2XzgWs_JR",
        "colab_type": "text"
      },
      "source": [
        "conv_autoencoder is from https://github.com/jellycsc/PyTorch-CIFAR-10-autoencoder/"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZpExFuC2VoX6",
        "colab_type": "text"
      },
      "source": [
        "同學們可以嘗試各種不同的VAE架構"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Wi8ds1fugCkR",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import torch\n",
        "from torch import nn\n",
        "import torch.nn.functional as F\n",
        "\n",
        "\n",
        "class fcn_autoencoder(nn.Module):\n",
        "    def __init__(self):\n",
        "        super(fcn_autoencoder, self).__init__()\n",
        "        self.encoder = nn.Sequential(\n",
        "            nn.Linear(32 * 32 * 3, 128),\n",
        "            nn.ReLU(True),\n",
        "            nn.Linear(128, 64),\n",
        "            nn.ReLU(True), nn.Linear(64, 12), nn.ReLU(True), nn.Linear(12, 3))\n",
        "        self.decoder = nn.Sequential(\n",
        "            nn.Linear(3, 12),\n",
        "            nn.ReLU(True),\n",
        "            nn.Linear(12, 64),\n",
        "            nn.ReLU(True),\n",
        "            nn.Linear(64, 128),\n",
        "            nn.ReLU(True), nn.Linear(128, 32 * 32 * 3\n",
        "            ), nn.Tanh())\n",
        "\n",
        "    def forward(self, x):\n",
        "        x = self.encoder(x)\n",
        "        x = self.decoder(x)\n",
        "        return x\n",
        "\n",
        "\n",
        "class conv_autoencoder(nn.Module):\n",
        "    def __init__(self):\n",
        "        super(conv_autoencoder, self).__init__()\n",
        "        self.encoder = nn.Sequential(\n",
        "            nn.Conv2d(3, 12, 4, stride=2, padding=1),            # [batch, 12, 16, 16]\n",
        "            nn.ReLU(),\n",
        "            nn.Conv2d(12, 24, 4, stride=2, padding=1),           # [batch, 24, 8, 8]\n",
        "            nn.ReLU(),\n",
        "\t\t\t      nn.Conv2d(24, 48, 4, stride=2, padding=1),           # [batch, 48, 4, 4]\n",
        "            nn.ReLU(),\n",
        "    # \t\t\tnn.Conv2d(48, 96, 4, stride=2, padding=1),           # [batch, 96, 2, 2]\n",
        "    #       nn.ReLU(),\n",
        "        )\n",
        "        self.decoder = nn.Sequential(\n",
        "#             nn.ConvTranspose2d(96, 48, 4, stride=2, padding=1),  # [batch, 48, 4, 4]\n",
        "#             nn.ReLU(),\n",
        "\t\t\t      nn.ConvTranspose2d(48, 24, 4, stride=2, padding=1),  # [batch, 24, 8, 8]\n",
        "            nn.ReLU(),\n",
        "\t\t\t      nn.ConvTranspose2d(24, 12, 4, stride=2, padding=1),  # [batch, 12, 16, 16]\n",
        "            nn.ReLU(),\n",
        "            nn.ConvTranspose2d(12, 3, 4, stride=2, padding=1),   # [batch, 3, 32, 32]\n",
        "            nn.Tanh(),\n",
        "        )\n",
        "\n",
        "    def forward(self, x):\n",
        "        x = self.encoder(x)\n",
        "        x = self.decoder(x)\n",
        "        return x\n",
        "\n",
        "\n",
        "class VAE(nn.Module):\n",
        "    def __init__(self):\n",
        "        super(VAE, self).__init__()\n",
        "\n",
        "        self.fc1 = nn.Linear(32*32*3, 400)\n",
        "        self.fc21 = nn.Linear(400, 20)\n",
        "        self.fc22 = nn.Linear(400, 20)\n",
        "        self.fc3 = nn.Linear(20, 400)\n",
        "        self.fc4 = nn.Linear(400, 32*32*3)\n",
        "\n",
        "    def encode(self, x):\n",
        "        h1 = F.relu(self.fc1(x))\n",
        "        return self.fc21(h1), self.fc22(h1)\n",
        "\n",
        "    def reparametrize(self, mu, logvar):\n",
        "        std = logvar.mul(0.5).exp_()\n",
        "        if torch.cuda.is_available():\n",
        "            eps = torch.cuda.FloatTensor(std.size()).normal_()\n",
        "        else:\n",
        "            eps = torch.FloatTensor(std.size()).normal_()\n",
        "        eps = Variable(eps)\n",
        "        return eps.mul(std).add_(mu)\n",
        "\n",
        "    def decode(self, z):\n",
        "        h3 = F.relu(self.fc3(z))\n",
        "        return F.sigmoid(self.fc4(h3))\n",
        "\n",
        "    def forward(self, x):\n",
        "        mu, logvar = self.encode(x)\n",
        "        z = self.reparametrize(mu, logvar)\n",
        "        return self.decode(z), mu, logvar\n",
        "\n",
        "\n",
        "def loss_vae(recon_x, x, mu, logvar, criterion):\n",
        "    \"\"\"\n",
        "    recon_x: generating images\n",
        "    x: origin images\n",
        "    mu: latent mean\n",
        "    logvar: latent log variance\n",
        "    \"\"\"\n",
        "    mse = criterion(recon_x, x)  # mse loss\n",
        "    # loss = 0.5 * sum(1 + log(sigma^2) - mu^2 - sigma^2)\n",
        "    KLD_element = mu.pow(2).add_(logvar.exp()).mul_(-1).add_(1).add_(logvar)\n",
        "    KLD = torch.sum(KLD_element).mul_(-0.5)\n",
        "    # KL divergence\n",
        "    return mse + KLD\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XKNUImqUhIeq",
        "colab_type": "text"
      },
      "source": [
        "# Training"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JoW1UrrxgI_U",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from torch.autograd import Variable\n",
        "from torch.utils.data import DataLoader\n",
        "from torch.optim import Adam, AdamW\n",
        "from torch.utils.data import (DataLoader, RandomSampler, SequentialSampler,\n",
        "                              TensorDataset)\n",
        "\n",
        "\n",
        "if task == 'ae':\n",
        "    num_epochs = 1000\n",
        "    batch_size = 128\n",
        "    learning_rate = 1e-3\n",
        "\n",
        "    #{'fcn', 'cnn', 'vae'} \n",
        "    model_type = 'cnn' \n",
        "\n",
        "    x = train\n",
        "    if model_type == 'fcn' or model_type == 'vae':\n",
        "        x = x.reshape(len(x), -1)\n",
        "        \n",
        "    data = torch.tensor(x, dtype=torch.float)\n",
        "    train_dataset = TensorDataset(data)\n",
        "    train_sampler = RandomSampler(train_dataset)\n",
        "    train_dataloader = DataLoader(train_dataset, sampler=train_sampler, batch_size=batch_size)\n",
        "\n",
        "\n",
        "    model_classes = {'fcn':fcn_autoencoder(), 'cnn':conv_autoencoder(), 'vae':VAE()}\n",
        "    model = model_classes[model_type].cuda()\n",
        "    criterion = nn.MSELoss()\n",
        "    optimizer = torch.optim.AdamW(\n",
        "        model.parameters(), lr=learning_rate)\n",
        "    \n",
        "    best_loss = np.inf\n",
        "    model.train()\n",
        "    for epoch in range(num_epochs):\n",
        "        for data in train_dataloader:\n",
        "            if model_type == 'cnn':\n",
        "                img = data[0].transpose(3, 1).cuda()\n",
        "            else:\n",
        "                img = data[0].cuda()\n",
        "            # ===================forward=====================\n",
        "            output = model(img)\n",
        "            if model_type == 'vae':\n",
        "                loss = loss_vae(output[0], img, output[1], output[2], criterion)\n",
        "            else:\n",
        "                loss = criterion(output, img)\n",
        "            # ===================backward====================\n",
        "            optimizer.zero_grad()\n",
        "            loss.backward()\n",
        "            optimizer.step()\n",
        "            # ===================save====================\n",
        "            if loss.item() < best_loss:\n",
        "                best_loss = loss.item()\n",
        "                torch.save(model, 'best_model_{}.pt'.format(model_type))\n",
        "        # ===================log========================\n",
        "        print('epoch [{}/{}], loss:{:.4f}'\n",
        "              .format(epoch + 1, num_epochs, loss.item()))\n",
        "        \n",
        "\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Wk0UxFuchLzR",
        "colab_type": "text"
      },
      "source": [
        "# Evaluation"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ku_wpuheWWNz",
        "colab_type": "text"
      },
      "source": [
        "將testing的圖片輸入model後，可以得到其重建的圖片，並對兩者取平方差。可以發現inlier的平方差應該與outlier的平方差形成差距明顯的兩群數據。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "r1PS_ApzhfOQ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "\n",
        "if task == 'ae':\n",
        "    if model_type == 'fcn' or model_type == 'vae':\n",
        "        y = test.reshape(len(test_tmp), -1)\n",
        "    else:\n",
        "        y = test\n",
        "        \n",
        "    data = torch.tensor(y, dtype=torch.float)\n",
        "    test_dataset = TensorDataset(data)\n",
        "    test_sampler = SequentialSampler(test_dataset)\n",
        "    test_dataloader = DataLoader(test_dataset, sampler=test_sampler, batch_size=batch_size)\n",
        "\n",
        "    model = torch.load('best_model_{}.pt'.format(model_type), map_location='cuda')\n",
        "\n",
        "    model.eval()\n",
        "    reconstructed = list()\n",
        "    for i, data in enumerate(test_dataloader): \n",
        "        if model_type == 'cnn':\n",
        "            img = data[0].transpose(3, 1).cuda()\n",
        "        else:\n",
        "            img = data[0].cuda()\n",
        "        output = model(img)\n",
        "        if model_type == 'cnn':\n",
        "            output = output.transpose(3, 1)\n",
        "        elif model_type == 'vae':\n",
        "            output = output[0]\n",
        "        reconstructed.append(output.cpu().detach().numpy())\n",
        "\n",
        "    reconstructed = np.concatenate(reconstructed, axis=0)\n",
        "    anomality = np.sqrt(np.sum(np.square(reconstructed - y).reshape(len(y), -1), axis=1))\n",
        "    y_pred = anomality\n",
        "    with open('prediction.csv', 'w') as f:\n",
        "        f.write('id,anomaly\\n')\n",
        "        for i in range(len(y_pred)):\n",
        "            f.write('{},{}\\n'.format(i+1, y_pred[i]))\n",
        "    # score = roc_auc_score(y_label, y_pred, average='micro')\n",
        "    # score = f1_score(y_label, y_pred, average='micro')\n",
        "    # print('auc score: {}'.format(score))\n"
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}